%% Script to load the .IBW files for measuring Vrest in Candellabrum Cells using Vramp.


% path='C:\Users\tomas\Documents\1. Harvard\2017 - Spring\Regehr lab rotation\Data\';
% path2 = '09_17\9_19_17\9_19_17_04\';
path=['C:\Users\Tom�s\Desktop\Temp\'];
D = IBWread([path 'RampAvg.ibw']);
trace=medfilt1(D.y, 5e2);
x=[1:length(trace)];
x2=[1:length(trace)];

figure; 
plot(x,trace);

index = (x > 2.5e3) & (x <6e3);
hold on;

shift= 2.5e3 ; %index is a logical vector,this helps figure out the starting point of the fitted part for plotting
plot(1+shift:sum(index)+shift,trace(index),'r','linewidth', 3) %plot in red the part of the trace that is specified by the variable index and will be used for linear fit

intersectRange= (x2>10e3) & (x<14.5e3); %part of the 
p=polyfit(x(index),trace(index)',1); %Make linear fit of the chosen indices
tracefit= p(1)*x +p(2);              %Use coefficients to make fit function
hold on;
plot(x,tracefit,'g','LineWidth', 1.2);               %plot the trace fit.
 
%Find index of the minimum between the fit and the trace (index of Vh in
%the function that describes the Voltage change (a linear function with
%negative slope.      

difference= abs(trace-tracefit');  % Absolute value of difference is important

intersectindex=find(difference == min(difference(intersectRange))); %Find the index of the intersection.
hold on;
plot(intersectindex,trace(intersectindex), 'k.','markers',12) %plot a black marker on the intersection.
ylim([(min(trace(1.5e4))-5) (max(trace(0.2e4))+5)]);
xlim([0.2e4 1.5e4]); %Ctrl + R this line to autoscale to full trace.
%axis tight;
xlabel('samples');
ylabel('pA');


%save(filename,variables)

voltCoefficients = polyfit ([1e3, 16e3], [100,-120],1);  % find linear fit coefficients for Vcommand
m2= voltCoefficients(1);
b2= voltCoefficients(2);

%hold on;
voltPlot=m2*x+b2;
%plot(x,voltPlot);
title([ path2(end-10) '.' path2(end-8:end-7) '.' path2(end-5:end-4) '  Cell ' path2(end-2:end-1) '  Vrest =' num2str(voltPlot(intersectindex)) 'mV']);
saveas(gcf,[path path2 path2(end-10) '.' path2(end-8:end-7) '.' path2(end-5:end-4) '  Cell ' path2(end-2:end-1)  'VRamp' ' Vrest =' num2str(voltPlot(intersectindex)) 'mV' '.eps'],'epsc');
saveas(gcf,[path path2 path2(end-10) '.' path2(end-8:end-7) '.' path2(end-5:end-4) '  Cell ' path2(end-2:end-1)  'VRamp' ' Vrest =' num2str(voltPlot(intersectindex)) 'mV' '.tif']);


 
%hold on 
%hax=axes
%SP=intersectindex; %your point goes here 
%line([SP SP],get(hax,'YLim'),'Color',[1 0 0])

